FUNCTION pythag_sp(a,b)
  
  USE nrtype
  IMPLICIT NONE
  
  REAL(SP), INTENT(IN) :: a,b
  REAL(SP) :: pythag_sp
  REAL(SP) :: absa,absb

  absa=abs(a)
  absb=abs(b)
  if (absa > absb) then
     pythag_sp=absa*sqrt(1.0_sp+(absb/absa)**2)
  else
     if (absb == 0.0) then
        pythag_sp=0.0
     else
        pythag_sp=absb*sqrt(1.0_sp+(absa/absb)**2)
     end if
  end if

END FUNCTION pythag_sp

FUNCTION pythag_dp(a,b)

  USE nrtype
  IMPLICIT NONE

  REAL(DP), INTENT(IN) :: a,b
  REAL(DP) :: pythag_dp
  REAL(DP) :: absa,absb

  absa=abs(a)
  absb=abs(b)
  if (absa > absb) then
     pythag_dp=absa*sqrt(1.0_dp+(absb/absa)**2)
  else
     if (absb == 0.0) then
        pythag_dp=0.0
     else
        pythag_dp=absb*sqrt(1.0_dp+(absa/absb)**2)
     end if
  end if

END FUNCTION pythag_dp
